Workshop Prep

Open https://github.com/dlab-geo/r-geospatial-workshop
- Download & unzip the zip file
- Make a note of the folder in which it is located
Start RStudio and open a new script file
Install required libraries in RStudio - if you do not have them already

install.packages(
  c("sp","rgdal","tmap","classInt","RColorBrewer",
    "ggplot2","leaflet", "ggmap"), dependencies=TRUE
)

Geospatial Data in R

Workshop Goals

Intro to working with geospatial data in R

geospatial data files and formats
Loading geospatial data in R
R packages for working with geospatial data
coordinate reference systems

Mapping geospatial data

Practice

About Me

Patty Frontiera

dlab.berkeley.edu

About you

Who are you?

Why are you here?

Follow Along

Open one of the tutorial files in a web browser

Slides r-geospatial-workshop-pt1.html

Tutorial r-geospatial-workshop-pt1-tutorial.html

Raw Code scripts/r-geospatial-workshop-pt1.R

Make sure you can cut and paste into RStudio

Geographic Data

are data about locations on or near the surface of the Earth.

Geospatial data

represent location more specifically with coordinates

46.130479, -117.134167

Coordinate Reference Systems

Coordinates only make sense when associated with a CRS!

Geographic Coordinates: Latitude and Longitude

Coordinate Reference Systems

Define:

the shape of the Earth
the origin (0,0 point)
the relationship between the system and the real world
the units

Because of variations in these, there are many geographic CRSs!

WGS84

The World Geodetic System of 1984 is the most widely used geographic coordinate reference system.

WGS84 is the default CRS for most GIS software

Almost all longitude and latitude data are assumed to be WGS84 unless otherwise specified

Historical data much trickier

Geospatial data are powerful!

You can

dynamically determine spatial metrics like area, length, distance and direction
spatial relationships like intersects, inside, contains, etc

and

link data by location, like census data and crime data

Spatial Data

Spatial data is a broader term than geographic data.

Methods for working with spatial data are valid for geospatial data

All spatial data are encoded with some type of coordinate reference system

Geospatial data require the CRS to be a model of locations on the Earth

Types of Spatial Data

Vector and Raster Data

Vector Data

Points, lines and Polygons

Raster Data

Regular grid of cells (or pixels)

We won’t be covering Raster data in this workshop*

Softare for working with Geospatial Data

Geospatial data require

software that can import, create, store, edit, visualize and analyze geospatial data

represented as geometric data objects referenced to the surface of the earth via CRSs
with methods to operate on those representations

GIS

We call software for working with geospatial data GIS

Geographic Information System

This term is commonly associated with desktop software applications.

Types of GIS Software

Desktop GIS - ArcGIS, QGIS

Spatial Databases - PostgreSQL/PostGIS

Web-based GIS - ArcGIS Online, CARTO

Software geospatial data support - Tableau

Programming languages with geospatial data support

R, Python, Javascript

Why R for Geospatial Data?

You already use R

Reproducibility

Free & Open Source

Strong support for geospatial data and analysis

Cutting edge

Geospatial Data in R

Prep reminder

Make sure you have the packages we are going to use installed and loaded

Make sure you have set your working directory to the location of the workshop files

install.packages(
  c("sp","rgdal","tmap","classInt","RColorBrewer",
    "ggplot2","leaflet", "ggmap"), dependencies=TRUE
)

Geospatial Data in R

There are many approaches to and packages for working with geospatial data in R.

One approach is to keep it simple and store geospatial data in a data frame.

This approach is most common when

the data are point data in CSV files and
you want to map rather than spatially transform or analyze the data

About the Sample Data

sf_properties_25ksample.csv

San Francisco Open Data Portal https://data.sfgov.org

SF Property Tax Rolls

This data set includes the Office of the Assessor-Recorder’s secured property tax roll spanning from 2007 to 2016.

We are using a subset of these data as a proxy for home values.

Load the CSV file into a data frame

sfhomes <- read.csv('data/sf_properties_25ksample.csv', 
                    stringsAsFactors = FALSE)

# Take a look at first 5 rows and a few of the columns
sfhomes[1:5,c("YearBuilt","totvalue","AreaSquareFeet","Neighborhood",
              "NumBedrooms")]

Make sure your working directory is set to the folder where you downloaded the workshop files!

sfhomes <- read.csv('data/sf_properties_25ksample.csv', 
                    stringsAsFactors = FALSE)

# Take a look at first 5 rows and a few of the columns
sfhomes[1:5,c("YearBuilt","totvalue","AreaSquareFeet","Neighborhood",
              "NumBedrooms")]

##   YearBuilt totvalue AreaSquareFeet       Neighborhood NumBedrooms
## 1      1923  1031391           2250 West of Twin Peaks           0
## 2      1909   775300           5830   Presidio Heights           3
## 3      1915  1116772           1350     Inner Richmond           2
## 4      1947   860130           1162    Sunset/Parkside           0
## 5      1981   322749           1463     Outer Richmond           3

Explore the data

class(sfhomes)            # what is the data object type?
dim(sfhomes)              # how many rows and columns
str(sfhomes)              # display the structure of the object
head(sfhomes)             # take a look at the first 10 records
summary(sfhomes)          # explore the range of values
summary(sfhomes$totvalue) # explore the range of values for one column
hist(sfhomes$totvalue)    # histogram for the totvalue column

Questions:

What columns contain the geographic data?
Are these data vector or raster data?
What type of geometry do the data contain?
- Points, lines, polygons, grid cells?
What is the CRS of these data?

Plot of points

Use the R base plot function to create a simple map

plot(sfhomes$lon, sfhomes$lat) # using base plot function

Plot of points

Use the R base plot function to create a simple map

plot(sfhomes$lon, sfhomes$lat) # using base plot function

`ggplot2`

Most widely used plotting library in R

Not specifically for geospatial data

But can be used to make fabulous maps

Great choice if you already know ggplot2

`ggplot2`

Load the library

library(ggplot2)

Maps with `ggplot2`

Basic map with ggplot

library(ggplot2)

ggplot() + geom_point(data=sfhomes, aes(lon,lat))

Maps with `ggplot2`

Basic map with ggplot

ggplot() + geom_point(data=sfhomes, aes(lon,lat), size=1)

Coord_map Option

ggplot() + geom_point(data=sfhomes, aes(lon,lat), size=1) + coord_map()

coord_map option

Allows you to associate a map projection with geographic coord data.

Map Projections

Map Projection: mathematial transformation from curved to flat surface

A Projected CRS applies a map projection to a Geographic CRS

Many Map Projections & Projected CRSs

All introduce distortion,

in shape, area, distance, direction, or combo
the larger the area the greater the distortion

No one map projection best for all purposes

Selection depends on location, extent and purpose

Different Projected CRSs

Map points symbolized by `totvalue`

Data driven symbology

ggplot() + geom_point(data=sfhomes, aes(lon,lat, col=totvalue)) + 
  coord_map()

Map points symbolized by `totvalue`

Data Order

What’s happening here?

sfhomes_low2high <- sfhomes[order(sfhomes$totvalue, decreasing = FALSE),]

ggplot() + 
  geom_point(data=sfhomes_low2high, aes(lon,lat, col=totvalue)) + 
  coord_map()

Try it - Does the output map look different from previous one?

Data Order

The order of the data in the data frame changes the map display!

Challenge

Map the sfhomes data in decreasing order by totvalue.

Decreasing order by `totvalue`

sfhomes_high2low <- sfhomes[order(sfhomes$totvalue, decreasing = T),]
ggplot() + geom_point(data=sfhomes_high2low, aes(lon,lat, col=totvalue)) + 
  coord_map()

More `ggplot` Goodness

What does this code do?

sfhomes2010_15 <- subset(sfhomes_low2high, as.numeric(SalesYear) > 2009)

ggplot() +
  geom_point(aes(lon, lat, col=totvalue), data = sfhomes2010_15 )  +
  facet_wrap(~ SalesYear)

More `ggplot` Goodness

Visual spatial analysis!

ggmap

`ggmap` extends `ggplot`

Create basemaps on which you can display your data.

Geocode place names and addresses to get point coordinates.

and more…

ggmap

Load the libary

library(ggmap)

## Google Maps API Terms of Service: http://developers.google.com/maps/terms.

## Please cite ggmap if you use it: see citation("ggmap") for details.

ggmap

Some ggmap functionality may require you to register a Google API key

Search - Google API Key
Create Google Maps API Key
Enable API key for geocoding, static maps API, javascript API

Register Google Maps Key

with ggmap

#register_google(key="AIzXXXXXXXXXXXXXXXXXxPSE") # your key here
register_google(key="AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo")

This key will be disabled after the workshop!

Re-install `ggmap`

If you have problems with ggmap you may need to reinstall the library from github.

See this StackOverflow discussion

#devtools::install_github("dkahle/ggmap")
#library(ggmap)

get_map

Fetch map data to plot using get_map

Use ?get_map for details

sf_map <- get_map("San Francisco, CA")

This use of the command leverages the Google Geocoding API to determine the coordinates of the named location.

get_map

Fetch map data to plot using get_map

Use ?get_map for details

sf_map <- get_map("San Francisco, CA")

## Source : https://maps.googleapis.com/maps/api/staticmap?center=San+Francisco,+CA&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=San%20Francisco%2C%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

Display the ggmap of SF

ggmap(sf_map)

Display the ggmap of SF

ggmap(sf_map)

`ggmap` with point overlay

Syntax similar to ggplot

# ggplot() +

ggmap(sf_map) + 
  geom_point(data=sfhomes, aes(x=lon, y=lat, col=totvalue))

`ggmap` with point overlay

ggmap(sf_map) +
  geom_point(data=sfhomes, aes(x=lon, y=lat, col=totvalue))

Customize `get_map`

Create a basemap zoomed to the extent of our data

# FIRST - subset the data
sfhomes15 <- subset(sfhomes, as.numeric(SalesYear) == 2015)

# Get the center point of the data
sf_ctr <- c(lon = mean(sfhomes$lon), lat = mean(sfhomes$lat))
sf_ctr  # take a look

# create the map
sf_basemap <- get_map(sf_ctr, zoom=12, scale=1)

Customize `get_map` Extent

Get a basemap zoomed to the extent of the 2015 data

# FIRST - subset the data
sfhomes15 <- subset(sfhomes, as.numeric(SalesYear) == 2015)

# Get the center point of the data
sf_ctr <- c(lon = mean(sfhomes15$lon), lat = mean(sfhomes15$lat))
sf_ctr  # take a look

##        lon        lat 
## -122.43289   37.75964

# create the map
sf_basemap <- get_map(sf_ctr, zoom=12, scale=1)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=37.759644,-122.432893&zoom=12&size=640x640&scale=1&maptype=terrain&language=en-EN&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

Challenge

Create a ggmap that uses the customized sf_basemap

Recreate our last map of property values (totvalue)

Map the points with `sf_basemap`

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat, col=totvalue))

Map Overlays

Finally, let’s add another geospatial data layer to our ggmap.

You can use this method to add as many layers as you want to a ggmap or ggplot.

Bart Stations and Landmarks

Use the read.csv function to read in a file of Bart Station locations.

bart <- read.csv("./data/bart.csv")
# take a look
head (bart)

##           X        Y             STATION OPERATOR DIST  CO
## 1 -122.2833 37.87406      NORTH BERKELEY     BART    4 ALA
## 2 -122.2682 37.86969   DOWNTOWN BERKELEY     BART    4 ALA
## 3 -122.2701 37.85321               ASHBY     BART    4 ALA
## 4 -122.2518 37.84451           ROCKRIDGE     BART    4 ALA
## 5 -122.2671 37.82871           MACARTHUR     BART    4 ALA
## 6 -122.2684 37.80808 19TH STREET/OAKLAND     BART    4 ALA

Add BART stations to map

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat))  +
  geom_point(data=bart, aes(x=X,y=Y), col="red")

Add BART stations to map

## Warning: Removed 34 rows containing missing values (geom_point).

Questons?

Geocoding with GGMAP

Geocoding

Determining the geographic coordinates for one or more addresses, zip codes, or named places.

Geocode SF Police Stations

police_stations <- read.csv("data/sf_police_addresses.csv", 
                            stringsAsFactors = F)
head(police_stations)

##   PoliceDistrict        Phone              Address          City State
## 1        Central 415-315-2400   766 Vallejo Street San Francisco    CA
## 2       Southern 415-575-6000      1251 3rd Street San Francisco    CA
## 3        Bayview 415-671-2300  201 Williams Avenue San Francisco    CA
## 4        Mission 415-558-5400  630 Valencia Street San Francisco    CA
## 5       Northern 415-614-3400 1125 Fillmore Street San Francisco    CA
## 6           Park 415-242-3000   1899 Waller Street San Francisco    CA

Geocode the Police Stations

#?geocode
geocode(police_stations)

Single column address

police_stations$full_address <- paste(police_stations$Address, 
                                      police_stations$City, police_stations$State)
head(police_stations)

##   PoliceDistrict        Phone              Address          City State
## 1        Central 415-315-2400   766 Vallejo Street San Francisco    CA
## 2       Southern 415-575-6000      1251 3rd Street San Francisco    CA
## 3        Bayview 415-671-2300  201 Williams Avenue San Francisco    CA
## 4        Mission 415-558-5400  630 Valencia Street San Francisco    CA
## 5       Northern 415-614-3400 1125 Fillmore Street San Francisco    CA
## 6           Park 415-242-3000   1899 Waller Street San Francisco    CA
##                            full_address
## 1   766 Vallejo Street San Francisco CA
## 2      1251 3rd Street San Francisco CA
## 3  201 Williams Avenue San Francisco CA
## 4  630 Valencia Street San Francisco CA
## 5 1125 Fillmore Street San Francisco CA
## 6   1899 Waller Street San Francisco CA

Geocode

station_coords <- geocode(police_stations$full_address)

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=766%20Vallejo%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1251%203rd%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=201%20Williams%20Avenue%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=630%20Valencia%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1125%20Fillmore%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1899%20Waller%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=461%206th%20Avenue%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1%20Sgt.%20John%20V.%20Young%20Lane%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=2345%2024th%20Avenue%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=301%20Eddy%20Street%20San%20Francisco%20CA&key=AIzaSyBE23PnI3bzQBJkhXuhXTenSfwwOVH6xUo

Add the lat/lons to the data

police_stations <- cbind(police_stations, station_coords)
police_stations[1:5,c("PoliceDistrict","full_address","lon","lat")]

##   PoliceDistrict                          full_address       lon      lat
## 1        Central   766 Vallejo Street San Francisco CA -122.4100 37.79872
## 2       Southern      1251 3rd Street San Francisco CA -122.3894 37.77238
## 3        Bayview  201 Williams Avenue San Francisco CA -122.3980 37.72976
## 4        Mission  630 Valencia Street San Francisco CA -122.4220 37.76285
## 5       Northern 1125 Fillmore Street San Francisco CA -122.4325 37.78022

Add the police stations

to your last map

Use a different color

Add Police Stations

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat))  +
  geom_point(data=bart, aes(x=X,y=Y), col="red", size=3) +
  geom_point(data=police_stations, aes(x=lon,y=lat), shape=22, 
             col="black", fill="grey", size=4)

Add Police Stations

## Warning: Removed 34 rows containing missing values (geom_point).

## About Google geocoding

You can geocode addresses or place names
Super accurate
You may need an API key
Limited to 2500 per day
Use ESRI Geocoding tools if you have more addresses

Let’s add one more layer

SF Landmarks

data/landmarks.csv

SF Landmarks

landmarks <- read.csv("./data/landmarks.csv")
head(landmarks)

##           X       Y id             name
## 1 -13626151 4551548  1       Coit Tower
## 2 -13630804 4544523  2       Twin Peaks
## 3 -13627654 4548289  3        City Hall
## 4 -13635101 4546904  4 Golden Gate Park

Map Landmarks

ggmap(sf_basemap) +
  geom_point(data=sfhomes15, aes(x=lon, y=lat))  +
  geom_point(data=bart, aes(x=X,y=Y), col="red", size=3) +
  geom_point(data=landmarks, aes(x=X,y=Y), shape=22, 
             col="black", fill="grey", size=4)

Map Landmarks

All good? If not, why not?

## Warning: Removed 34 rows containing missing values (geom_point).

## Warning: Removed 4 rows containing missing values (geom_point).

GGMap and GGPlot are great!

BUT!

There are limits to what you can do with geospatial data stored in a data frame.

Can’t transform Coordinate Data

These data are not in geographic coordinates - longitude and latitude.

You can’t map these with ggmap.

##           X       Y id             name
## 1 -13626151 4551548  1       Coit Tower
## 2 -13630804 4544523  2       Twin Peaks
## 3 -13627654 4548289  3        City Hall
## 4 -13635101 4546904  4 Golden Gate Park

Can’t read & plot geospatial data files

The ESRI Shapefile is the most common file format for geospatial data.

Can’t read in or directly plot data in a shapefile with ggplot or ggmap.

Can’t Perform Spatial Analysis

Spatial data objects and methods are needed to answer questions like

What properties are in the Noe Valley neighborhood?
What is the average property value in each SF neighborhood?
What is the area of each SF Neighborhood and the property density?
What properties are within walking distance (.25 miles) of the Mission neighborhood?

Spatial Data Objects in R

sp Package

The `sp` Package

Classes and Methods for Spatial Data

The SP package is most commonly used to construct and manipulate spatial data objects in R.

Hundreds of other R packages that do things with spatial data typically build on SP objects.

sp package

`sp` package

Load the library sp

Take a look at the different types of spatial object classes supported by sp

library(sp)
getClass("Spatial")

sp

## Warning: package 'sp' was built under R version 3.4.3

## Class "Spatial" [package "sp"]
## 
## Slots:
##                               
## Name:         bbox proj4string
## Class:      matrix         CRS
## 
## Known Subclasses: 
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2

sp Vector Objects

Geometry	Spatial Object	Spatial Object with Attributes
Points	SpatialPoints	SpatialPointsDataFrame
Lines	SpatialLines	SpatialLinesDataFrame
Polygons	SpatialPolygons	SpatialPolygonsDataFrame

We use the S*DF objects most frequently!

From Data Frame to SpatialPointsDataFrame

Let’s transform the sfhomes15 data frame to an sp object of type SpatialPointsDataFrame

`sp::coordinates()`

Use the sp::coordinates() method

Sets or retrieves spatial coordinates

When transforming a DF to SPDF you are setting the coordinates.

Note - the transformation happens in place!

Data Frame to SPDF

sfhomes15_sp <- sfhomes15  # Make a copy - why?

class(sfhomes15_sp)  # check the class of the object

coordinates(sfhomes15_sp) <- c('lon','lat')  #  Make it spatial - ORDER MATTERS!!

class(sfhomes15_sp)  # check the class of the object

Data Frame to SPDF

sfhomes15_sp <- sfhomes15  # Make a copy - why?

class(sfhomes15_sp)  # check the class of the object

## [1] "data.frame"

coordinates(sfhomes15_sp) <- c('lon','lat')  #  Make it spatial - ORDER MATTERS!!

class(sfhomes15_sp)  # check the class of the object

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

sp::coordinates()

You transformed a data frame to an SPDF using:

coordinates(sfhomes15_sp) <- c('lon','lat')

Now try this:

coordinates(sfhomes15_sp)

sp::coordinates()

coordinates(sfhomes15_sp)

##             lon      lat
## 24    -122.4741 37.73601
## 35    -122.3965 37.75920
## 76    -122.3942 37.77181
## 79    -122.4754 37.71518
## 85    -122.4210 37.75876
## 90    -122.4149 37.79362
## 92    -122.4370 37.74020
## 98    -122.4437 37.78881
## 223   -122.4754 37.71492
## 229   -122.3919 37.78868
## 273   -122.3968 37.77945
## 341   -122.5040 37.75515
## 363   -122.4041 37.76422
## 366   -122.3934 37.77764
## 388   -122.4166 37.74931
## 451   -122.4112 37.73683
## 467   -122.4507 37.76448
## 473   -122.3990 37.75450
## 487   -122.4195 37.76649
## 534   -122.4225 37.78676
## 548   -122.4410 37.70854
## 549   -122.5009 37.74485
## 550   -122.4220 37.73048
## 559   -122.4947 37.75832
## 567   -122.3942 37.77181
## 595   -122.4187 37.80469
## 624   -122.4019 37.77323
## 631   -122.4229 37.80586
## 646   -122.4225 37.78676
## 658   -122.3919 37.78868
## 665   -122.4261 37.74349
## 668   -122.3904 37.78139
## 683   -122.4458 37.74398
## 774   -122.3883 37.73754
## 776   -122.4767 37.75347
## 864   -122.4326 37.74229
## 974   -122.4380 37.72624
## 1018  -122.4452 37.76906
## 1097  -122.4937 37.77999
## 1155  -122.5047 37.75390
## 1164  -122.4498 37.75451
## 1167  -122.4055 37.80424
## 1188  -122.4194 37.75596
## 1242  -122.4088 37.77953
## 1273  -122.4331 37.76019
## 1374  -122.3918 37.77899
## 1379  -122.4070 37.75674
## 1487  -122.4236 37.77800
## 1598  -122.4758 37.75497
## 1600  -122.4419 37.74582
## 1618  -122.5030 37.74587
## 1668  -122.3989 37.78990
## 1698  -122.4236 37.77800
## 1747  -122.3972 37.76314
## 1755  -122.4434 37.75159
## 1768  -122.4242 37.80551
## 1804  -122.4549 37.76402
## 1828  -122.4320 37.78195
## 1833  -122.4405 37.71825
## 1913  -122.4146 37.74714
## 1999  -122.3942 37.77181
## 2024  -122.4676 37.71048
## 2033  -122.4249 37.76548
## 2066  -122.4243 37.75206
## 2073  -122.4097 37.77613
## 2096  -122.3996 37.73859
## 2248  -122.4045 37.79929
## 2258  -122.4442 37.78322
## 2285  -122.3919 37.78868
## 2314  -122.4026 37.77879
## 2329  -122.4368 37.78228
## 2339  -122.4240 37.74366
## 2355  -122.4126 37.79035
## 2416  -122.4579 37.72632
## 2428  -122.4089 37.79286
## 2436  -122.4409 37.80451
## 2443  -122.4127 37.78995
## 2500  -122.4357 37.76468
## 2508  -122.4666 37.76347
## 2537  -122.4870 37.74605
## 2538  -122.3942 37.77181
## 2627  -122.4219 37.77389
## 2692  -122.3958 37.77705
## 2717  -122.4116 37.77806
## 2734  -122.4290 37.72977
## 2793  -122.4764 37.73390
## 2914  -122.4183 37.76712
## 2969  -122.4711 37.72241
## 2975  -122.4393 37.78018
## 3127  -122.4393 37.78951
## 3157  -122.3702 37.72838
## 3205  -122.4398 37.78746
## 3262  -122.4540 37.78251
## 3288  -122.4698 37.75298
## 3296  -122.4607 37.73104
## 3306  -122.4081 37.77962
## 3354  -122.4849 37.74663
## 3403  -122.4209 37.75792
## 3436  -122.4854 37.74788
## 3471  -122.4856 37.78116
## 3485  -122.4177 37.80345
## 3516  -122.4354 37.72201
## 3578  -122.4199 37.77606
## 3611  -122.4269 37.80229
## 3658  -122.4108 37.76534
## 3754  -122.3888 37.72872
## 3772  -122.3891 37.73264
## 4023  -122.4827 37.77762
## 4025  -122.4450 37.76430
## 4058  -122.4328 37.71546
## 4076  -122.4660 37.74828
## 4111  -122.5023 37.74333
## 4114  -122.4970 37.74468
## 4140  -122.4741 37.73601
## 4150  -122.4060 37.79853
## 4249  -122.4302 37.74726
## 4270  -122.4194 37.76710
## 4287  -122.4676 37.71048
## 4369  -122.4950 37.73327
## 4376  -122.3919 37.78868
## 4393  -122.3972 37.77570
## 4427  -122.4332 37.77869
## 4466  -122.4786 37.77394
## 4469  -122.4210 37.76470
## 4483  -122.4318 37.71761
## 4514  -122.4214 37.78139
## 4526  -122.4153 37.74047
## 4527  -122.4209 37.79368
## 4555  -122.4145 37.79167
## 4561  -122.3910 37.78938
## 4566  -122.4031 37.79058
## 4578  -122.3906 37.75501
## 4579  -122.4154 37.74537
## 4592  -122.3922 37.78607
## 4599  -122.4385 37.75980
## 4602  -122.4785 37.76295
## 4618  -122.4137 37.79099
## 4648  -122.4635 37.76461
## 4663  -122.4160 37.74563
## 4669  -122.3918 37.77899
## 4759  -122.4911 37.75318
## 4800  -122.4137 37.79099
## 4825  -122.4147 37.72528
## 4853  -122.4388 37.75978
## 4864  -122.3696 37.72900
## 4896  -122.4245 37.79493
## 4899  -122.3902 37.75236
## 4921  -122.5103 37.77369
## 4939  -122.3930 37.73299
## 4944  -122.3919 37.78868
## 4970  -122.4646 37.77383
## 4986  -122.3904 37.78139
## 5053  -122.4287 37.77484
## 5085  -122.4244 37.78587
## 5104  -122.4255 37.72744
## 5108  -122.4364 37.77781
## 5109  -122.4394 37.80326
## 5114  -122.4853 37.75287
## 5188  -122.4302 37.74806
## 5228  -122.3881 37.73174
## 5247  -122.4242 37.80551
## 5252  -122.4296 37.79988
## 5261  -122.4048 37.79879
## 5287  -122.4335 37.79946
## 5288  -122.4611 37.73481
## 5316  -122.3918 37.77899
## 5320  -122.4481 37.77945
## 5321  -122.4477 37.71766
## 5334  -122.3919 37.78868
## 5350  -122.4336 37.72535
## 5351  -122.4540 37.78251
## 5360  -122.4264 37.72624
## 5366  -122.3922 37.78607
## 5398  -122.3905 37.78759
## 5425  -122.3919 37.78868
## 5430  -122.4263 37.72116
## 5450  -122.4175 37.79865
## 5461  -122.3989 37.78577
## 5467  -122.4175 37.77341
## 5511  -122.3919 37.78868
## 5519  -122.4093 37.74024
## 5521  -122.4816 37.77294
## 5522  -122.4194 37.75596
## 5547  -122.3959 37.79054
## 5561  -122.4213 37.78884
## 5571  -122.4194 37.76710
## 5609  -122.4245 37.79920
## 5647  -122.4085 37.75911
## 5654  -122.4673 37.75439
## 5671  -122.4459 37.75811
## 5691  -122.4806 37.77724
## 5705  -122.4167 37.76091
## 5712  -122.3702 37.72838
## 5732  -122.4207 37.73596
## 5740  -122.4101 37.77808
## 5757  -122.4752 37.71513
## 5765  -122.3922 37.73203
## 5827  -122.4036 37.71724
## 5978  -122.4264 37.73053
## 5997  -122.3886 37.77159
## 6043  -122.5068 37.77785
## 6056  -122.4089 37.79286
## 6134  -122.4191 37.70898
## 6140  -122.4107 37.76483
## 6181  -122.3974 37.78575
## 6191  -122.4321 37.73585
## 6208  -122.4774 37.75583
## 6252  -122.4355 37.74743
## 6258  -122.3904 37.78271
## 6342  -122.4846 37.75914
## 6365  -122.4356 37.75016
## 6410  -122.3824 37.72887
## 6430  -122.4699 37.75813
## 6477  -122.4656 37.71671
## 6581  -122.4792 37.77343
## 6687  -122.4243 37.79352
## 6691  -122.4371 37.71173
## 6694  -122.4399 37.71837
## 6721  -122.4308 37.73647
## 6723  -122.4566 37.75242
## 6782  -122.4298 37.71104
## 6787  -122.4405 37.80504
## 6814  -122.4848 37.73940
## 6863  -122.4122 37.79293
## 6892  -122.4419 37.74582
## 6903  -122.4084 37.79154
## 6915  -122.4144 37.77516
## 6918  -122.3958 37.77816
## 6932  -122.4244 37.74805
## 6933  -122.4368 37.78249
## 6988  -122.5064 37.77555
## 7006  -122.4528 37.76708
## 7009  -122.4947 37.76067
## 7032  -122.4050 37.77762
## 7073  -122.4194 37.75596
## 7142  -122.4953 37.73096
## 7149  -122.4104 37.75357
## 7175  -122.4370 37.76934
## 7206  -122.4936 37.77938
## 7303  -122.4327 37.77299
## 7320  -122.4194 37.76710
## 7337  -122.4752 37.71551
## 7395  -122.4263 37.76856
## 7404  -122.3939 37.75915
## 7447  -122.4113 37.72850
## 7462  -122.4127 37.76517
## 7493  -122.4194 37.75596
## 7533  -122.4515 37.75418
## 7555  -122.4071 37.77553
## 7564  -122.4389 37.73129
## 7600  -122.4097 37.77613
## 7649  -122.4290 37.78872
## 7651  -122.3958 37.77816
## 7704  -122.4337 37.76214
## 7735  -122.3904 37.78271
## 7832  -122.4853 37.75355
## 7833  -122.3983 37.75753
## 7848  -122.4490 37.76795
## 7936  -122.4087 37.77310
## 7981  -122.4489 37.71967
## 7991  -122.4325 37.77805
## 8047  -122.4019 37.78005
## 8089  -122.4194 37.75596
## 8127  -122.4685 37.76188
## 8161  -122.4307 37.74163
## 8275  -122.4126 37.77870
## 8322  -122.3924 37.73700
## 8344  -122.4126 37.77870
## 8378  -122.4900 37.78324
## 8401  -122.4642 37.71638
## 8410  -122.4279 37.77408
## 8419  -122.4016 37.75145
## 8422  -122.4836 37.71006
## 8438  -122.4668 37.73426
## 8517  -122.3919 37.78868
## 8537  -122.3980 37.71952
## 8590  -122.4291 37.73289
## 8612  -122.4194 37.75596
## 8652  -122.4222 37.72570
## 8694  -122.3919 37.78868
## 8703  -122.3958 37.77339
## 8713  -122.4301 37.77286
## 8730  -122.4126 37.77870
## 8818  -122.3897 37.78278
## 8823  -122.4361 37.78937
## 8842  -122.4257 37.79707
## 8844  -122.4259 37.78673
## 8888  -122.4318 37.79685
## 8891  -122.4277 37.77402
## 8899  -122.3896 37.78063
## 8907  -122.4728 37.75198
## 8932  -122.4102 37.74514
## 8962  -122.3942 37.77181
## 8983  -122.4206 37.80236
## 9049  -122.4651 37.75119
## 9054  -122.4371 37.73361
## 9179  -122.4262 37.72086
## 9216  -122.4160 37.73419
## 9219  -122.4550 37.74285
## 9356  -122.4522 37.73356
## 9376  -122.3982 37.75309
## 9384  -122.4152 37.73405
## 9469  -122.4174 37.80442
## 9479  -122.4251 37.72966
## 9490  -122.4835 37.74208
## 9569  -122.4723 37.73269
## 9639  -122.4883 37.77411
## 9640  -122.4509 37.73126
## 9647  -122.3909 37.75967
## 9675  -122.3891 37.76366
## 9691  -122.4381 37.78210
## 9764  -122.4875 37.76255
## 9776  -122.4422 37.72376
## 9806  -122.3962 37.77393
## 9830  -122.4642 37.71673
## 9906  -122.4433 37.76135
## 9959  -122.3934 37.77764
## 9966  -122.4639 37.72971
## 9985  -122.4411 37.75920
## 10010 -122.4279 37.76639
## 10046 -122.4755 37.71495
## 10063 -122.3879 37.73919
## 10066 -122.4429 37.78894
## 10100 -122.3890 37.73260
## 10157 -122.4377 37.79704
## 10164 -122.4180 37.74667
## 10170 -122.4453 37.79137
## 10229 -122.4137 37.79099
## 10239 -122.4451 37.78484
## 10257 -122.4102 37.74590
## 10263 -122.4353 37.75147
## 10272 -122.4974 37.78118
## 10273 -122.4401 37.75541
## 10307 -122.4167 37.79138
## 10338 -122.4343 37.73743
## 10343 -122.4108 37.77597
## 10345 -122.5037 37.74846
## 10351 -122.3910 37.78938
## 10353 -122.4865 37.76343
## 10410 -122.3940 37.75300
## 10415 -122.4287 37.73629
## 10419 -122.4344 37.76357
## 10499 -122.4531 37.73537
## 10530 -122.4150 37.73959
## 10578 -122.4594 37.73954
## 10583 -122.4760 37.71558
## 10600 -122.4832 37.78375
## 10607 -122.4911 37.73160
## 10632 -122.4347 37.78713
## 10646 -122.4518 37.76347
## 10661 -122.4066 37.73036
## 10781 -122.4124 37.80434
## 10782 -122.4554 37.71773
## 10800 -122.4152 37.73705
## 10857 -122.4178 37.76399
## 10887 -122.4260 37.74559
## 10914 -122.4349 37.74579
## 10957 -122.4173 37.74421
## 10967 -122.4587 37.73858
## 10985 -122.4325 37.77540
## 11005 -122.4558 37.73265
## 11023 -122.4173 37.73753
## 11035 -122.4386 37.76810
## 11094 -122.4275 37.72671
## 11158 -122.3937 37.78702
## 11160 -122.4236 37.77800
## 11170 -122.5054 37.76363
## 11171 -122.4001 37.79803
## 11180 -122.4757 37.71556
## 11181 -122.4359 37.75531
## 11202 -122.4133 37.78070
## 11258 -122.4437 37.80053
## 11265 -122.4865 37.76465
## 11284 -122.4117 37.75704
## 11286 -122.4194 37.75596
## 11305 -122.4136 37.79298
## 11385 -122.4424 37.74562
## 11451 -122.4370 37.73761
## 11498 -122.4438 37.71852
## 11554 -122.4990 37.77913
## 11564 -122.4194 37.75596
## 11576 -122.3910 37.78938
## 11598 -122.3824 37.72852
## 11603 -122.4748 37.71490
## 11622 -122.3988 37.71372
## 11636 -122.4678 37.72225
## 11690 -122.4415 37.78044
## 11728 -122.4793 37.75869
## 11729 -122.4165 37.73373
## 11808 -122.4457 37.75769
## 11830 -122.4965 37.77600
## 11844 -122.4468 37.73736
## 11873 -122.4474 37.73477
## 11876 -122.4126 37.79035
## 11882 -122.3970 37.75975
## 11918 -122.4438 37.76155
## 11927 -122.4333 37.73432
## 11968 -122.3936 37.78370
## 11969 -122.3907 37.78117
## 11971 -122.4742 37.75636
## 11972 -122.4542 37.77443
## 12052 -122.4385 37.75442
## 12232 -122.5060 37.73665
## 12285 -122.4786 37.77955
## 12350 -122.4672 37.73270
## 12383 -122.4093 37.77680
## 12428 -122.4274 37.77132
## 12467 -122.3971 37.78091
## 12533 -122.4024 37.77546
## 12559 -122.4130 37.80397
## 12569 -122.4114 37.72936
## 12622 -122.5059 37.73682
## 12668 -122.4131 37.79777
## 12672 -122.4154 37.73877
## 12675 -122.4357 37.73740
## 12730 -122.4011 37.73263
## 12770 -122.4221 37.72534
## 12783 -122.3910 37.78938
## 12865 -122.3919 37.78868
## 12866 -122.4891 37.75507
## 13000 -122.3968 37.77945
## 13015 -122.4341 37.76397
## 13054 -122.4165 37.77528
## 13177 -122.4043 37.78622
## 13251 -122.4294 37.77067
## 13252 -122.3933 37.77891
## 13268 -122.4316 37.76576
## 13293 -122.4936 37.77938
## 13324 -122.4116 37.75238
## 13382 -122.4513 37.76809
## 13391 -122.4248 37.80057
## 13414 -122.4433 37.78314
## 13422 -122.4141 37.80393
## 13470 -122.4439 37.71527
## 13472 -122.4815 37.74502
## 13478 -122.4214 37.78139
## 13507 -122.4298 37.73805
## 13529 -122.4242 37.80551
## 13550 -122.4439 37.80384
## 13581 -122.3956 37.73468
## 13593 -122.4286 37.74623
## 13595 -122.4617 37.78577
## 13662 -122.3918 37.77899
## 13675 -122.3696 37.72900
## 13701 -122.3962 37.77393
## 13722 -122.4025 37.73527
## 13777 -122.4750 37.71499
## 13821 -122.4032 37.78824
## 13895 -122.4485 37.73998
## 13903 -122.3919 37.78868
## 13926 -122.4141 37.75915
## 13936 -122.4296 37.73278
## 14044 -122.4484 37.76533
## 14072 -122.4546 37.74500
## 14087 -122.4642 37.75929
## 14239 -122.4665 37.71392
## 14241 -122.4456 37.80221
## 14256 -122.4990 37.77356
## 14287 -122.4234 37.77253
## 14338 -122.4475 37.77957
## 14345 -122.4192 37.75400
## 14353 -122.3696 37.72900
## 14367 -122.4253 37.79247
## 14413 -122.4397 37.76158
## 14496 -122.4276 37.79552
## 14517 -122.4617 37.75037
## 14679 -122.4565 37.72132
## 14690 -122.4365 37.71167
## 14699 -122.4375 37.76305
## 14721 -122.3919 37.78868
## 14728 -122.3934 37.77764
## 14769 -122.4258 37.78636
## 14773 -122.4326 37.75640
## 14777 -122.3936 37.78370
## 14803 -122.4121 37.77255
## 14870 -122.4481 37.73958
## 14965 -122.4093 37.74690
## 14974 -122.5070 37.78027
## 15069 -122.4340 37.74511
## 15076 -122.4733 37.77618
## 15143 -122.4234 37.77253
## 15159 -122.4937 37.78119
## 15167 -122.4129 37.78791
## 15179 -122.4668 37.74374
## 15186 -122.4476 37.77794
## 15187 -122.3991 37.78318
## 15197 -122.3943 37.77898
## 15208 -122.4557 37.71642
## 15217 -122.4678 37.72282
## 15260 -122.4416 37.71429
## 15267 -122.4507 37.74389
## 15319 -122.4345 37.78565
## 15327 -122.4957 37.77863
## 15328 -122.4564 37.70934
## 15349 -122.4151 37.73756
## 15357 -122.3761 37.73193
## 15426 -122.4009 37.76461
## 15474 -122.4009 37.76461
## 15485 -122.5046 37.75155
## 15486 -122.4718 37.76104
## 15498 -122.4174 37.74938
## 15525 -122.4124 37.77093
## 15551 -122.3943 37.78613
## 15567 -122.4055 37.80424
## 15580 -122.3914 37.78451
## 15587 -122.5046 37.74607
## 15589 -122.4339 37.72119
## 15604 -122.4248 37.75359
## 15615 -122.4098 37.71496
## 15621 -122.4667 37.73841
## 15666 -122.4379 37.75930
## 15669 -122.4194 37.75596
## 15713 -122.3702 37.72838
## 15717 -122.4382 37.75646
## 15790 -122.4021 37.76545
## 15887 -122.4311 37.76352
## 15918 -122.3942 37.77181
## 15933 -122.4178 37.76399
## 15945 -122.4390 37.79617
## 16033 -122.4014 37.78629
## 16076 -122.4149 37.79362
## 16105 -122.4489 37.76179
## 16173 -122.4135 37.72849
## 16199 -122.4014 37.77974
## 16204 -122.4327 37.76614
## 16210 -122.4210 37.75876
## 16212 -122.4126 37.77870
## 16215 -122.3958 37.77705
## 16229 -122.4400 37.71831
## 16310 -122.3988 37.72731
## 16315 -122.4254 37.80339
## 16389 -122.4800 37.75339
## 16440 -122.4970 37.74495
## 16450 -122.4194 37.75596
## 16458 -122.4624 37.77409
## 16474 -122.4542 37.77443
## 16480 -122.4014 37.77974
## 16500 -122.4126 37.77870
## 16507 -122.4300 37.78804
## 16517 -122.4240 37.76628
## 16531 -122.4266 37.73311
## 16549 -122.5038 37.74097
## 16640 -122.4384 37.77517
## 16643 -122.4517 37.74504
## 16712 -122.4194 37.75596
## 16755 -122.3895 37.73005
## 16778 -122.4708 37.75810
## 16796 -122.4035 37.72002
## 16800 -122.3696 37.72900
## 16802 -122.4042 37.71115
## 16921 -122.3696 37.72900
## 16934 -122.4054 37.77794
## 17006 -122.4459 37.77465
## 17051 -122.4560 37.75650
## 17073 -122.4386 37.73261
## 17084 -122.3942 37.77181
## 17091 -122.4750 37.71499
## 17156 -122.4286 37.71810
## 17168 -122.4766 37.76075
## 17199 -122.3910 37.78938
## 17200 -122.4151 37.78524
## 17223 -122.4011 37.77047
## 17249 -122.4159 37.74769
## 17252 -122.4908 37.73143
## 17265 -122.4043 37.80443
## 17304 -122.4014 37.77974
## 17345 -122.4245 37.79493
## 17358 -122.4704 37.72549
## 17384 -122.3904 37.78271
## 17404 -122.4080 37.77576
## 17410 -122.3904 37.75368
## 17465 -122.5033 37.74841
## 17491 -122.4015 37.73255
## 17510 -122.4359 37.75497
## 17517 -122.4519 37.76025
## 17542 -122.4124 37.77093
## 17562 -122.4514 37.70881
## 17613 -122.3893 37.73988
## 17618 -122.4249 37.75375
## 17698 -122.4262 37.79417
## 17715 -122.4140 37.75590
## 17754 -122.4360 37.80400
## 17767 -122.4916 37.78122
## 17857 -122.4065 37.77812
## 17867 -122.4777 37.76001
## 17876 -122.5107 37.77399
## 17881 -122.4226 37.73575
## 17887 -122.4188 37.77159
## 17912 -122.3886 37.77159
## 17929 -122.4289 37.72161
## 17933 -122.4452 37.71900
## 17936 -122.4458 37.72772
## 17986 -122.4429 37.73288
## 17988 -122.4473 37.73331
## 18002 -122.4422 37.77803
## 18018 -122.4772 37.74409
## 18050 -122.4198 37.76161
## 18053 -122.4038 37.76405
## 18059 -122.4221 37.74746
## 18071 -122.4271 37.73809
## 18107 -122.4388 37.73398
## 18164 -122.4160 37.73504
## 18227 -122.4805 37.75531
## 18231 -122.4671 37.76370
## 18323 -122.4229 37.80586
## 18363 -122.3727 37.73057
## 18386 -122.4767 37.77714
## 18433 -122.4408 37.78769
## 18439 -122.4345 37.73587
## 18444 -122.4716 37.78657
## 18453 -122.3928 37.78205
## 18470 -122.4682 37.75544
## 18504 -122.5075 37.74854
## 18511 -122.3959 37.79054
## 18553 -122.3942 37.77181
## 18584 -122.4140 37.75590
## 18619 -122.4315 37.73498
## 18627 -122.4716 37.75516
## 18642 -122.4780 37.78258
## 18653 -122.4258 37.72884
## 18718 -122.4391 37.75547
## 18733 -122.4294 37.72902
## 18742 -122.4194 37.75596
## 18792 -122.4181 37.78824
## 18795 -122.3938 37.78638
## 18801 -122.4288 37.72584
## 18871 -122.4234 37.77253
## 18886 -122.4305 37.74402
## 18913 -122.4280 37.74869
## 18960 -122.4387 37.71054
## 18968 -122.4041 37.80347
## 18971 -122.4058 37.80159
## 18990 -122.4398 37.75879
## 19015 -122.4613 37.76457
## 19065 -122.4121 37.77187
## 19082 -122.4684 37.76042
## 19105 -122.4836 37.71006
## 19143 -122.4642 37.78458
## 19153 -122.4520 37.73468
## 19197 -122.4102 37.75679
## 19202 -122.4412 37.71904
## 19227 -122.3896 37.78063
## 19281 -122.4309 37.72308
## 19325 -122.4287 37.77642
## 19375 -122.4152 37.73366
## 19532 -122.4596 37.78655
## 19535 -122.5019 37.77481
## 19570 -122.4623 37.73544
## 19626 -122.4348 37.74426
## 19641 -122.3919 37.78868
## 19718 -122.4194 37.75596
## 19765 -122.4194 37.75596
## 19782 -122.3911 37.73912
## 19802 -122.4152 37.79845
## 19811 -122.3951 37.72508
## 19832 -122.4798 37.78437
## 19838 -122.4301 37.74120
## 19845 -122.4435 37.75964
## 19880 -122.4243 37.79352
## 19897 -122.3934 37.77764
## 19932 -122.4590 37.71972
## 19933 -122.4261 37.73585
## 19954 -122.4934 37.76073
## 19991 -122.4242 37.79788
## 20009 -122.4528 37.71555
## 20026 -122.3909 37.75967
## 20061 -122.4382 37.80320
## 20095 -122.3919 37.78868
## 20102 -122.4936 37.74282
## 20221 -122.4251 37.75560
## 20235 -122.4709 37.75704
## 20293 -122.4937 37.77999
## 20310 -122.4252 37.74835
## 20320 -122.4355 37.71758
## 20321 -122.4739 37.74485
## 20347 -122.5013 37.76014
## 20356 -122.4517 37.73356
## 20426 -122.4014 37.77974
## 20484 -122.4345 37.80458
## 20562 -122.3974 37.78575
## 20605 -122.4729 37.77471
## 20611 -122.3922 37.78607
## 20655 -122.4621 37.77642
## 20687 -122.4548 37.78306
## 20743 -122.4469 37.71901
## 20754 -122.4493 37.77461
## 20768 -122.3962 37.77393
## 20778 -122.4743 37.76473
## 20779 -122.4024 37.76117
## 20804 -122.5022 37.74887
## 20871 -122.3922 37.73197
## 20877 -122.4180 37.76654
## 20928 -122.4731 37.77999
## 20971 -122.4816 37.77795
## 20979 -122.4582 37.73860
## 21007 -122.3696 37.72900
## 21053 -122.4194 37.75596
## 21064 -122.4263 37.76856
## 21079 -122.4425 37.80069
## 21093 -122.4336 37.74361
## 21105 -122.4245 37.79920
## 21124 -122.4399 37.79230
## 21152 -122.4642 37.72043
## 21174 -122.4227 37.79931
## 21201 -122.4340 37.77121
## 21231 -122.4272 37.72286
## 21255 -122.4094 37.75902
## 21374 -122.4208 37.73601
## 21424 -122.3702 37.72838
## 21455 -122.4195 37.80612
## 21456 -122.3932 37.78006
## 21458 -122.4980 37.75912
## 21479 -122.4331 37.79233
## 21486 -122.3922 37.78607
## 21496 -122.4739 37.78414
## 21499 -122.4031 37.79058
## 21508 -122.4606 37.72592
## 21522 -122.4188 37.77159
## 21599 -122.4630 37.74288
## 21626 -122.4271 37.78448
## 21632 -122.4368 37.74525
## 21644 -122.4176 37.76271
## 21652 -122.3901 37.75370
## 21678 -122.3702 37.72838
## 21694 -122.4414 37.79054
## 21711 -122.4126 37.77870
## 21712 -122.3904 37.78139
## 21790 -122.4446 37.79053
## 21791 -122.5024 37.73799
## 21834 -122.4218 37.79478
## 21889 -122.4157 37.79792
## 21932 -122.4520 37.71732
## 21957 -122.4282 37.79869
## 21972 -122.4552 37.74086
## 21985 -122.4549 37.73345
## 22007 -122.4195 37.74526
## 22024 -122.5006 37.75075
## 22064 -122.4569 37.75564
## 22085 -122.4497 37.72884
## 22111 -122.3933 37.77891
## 22167 -122.5041 37.73937
## 22254 -122.3992 37.79825
## 22255 -122.4107 37.77166
## 22273 -122.4448 37.71416
## 22324 -122.3986 37.71395
## 22340 -122.4376 37.75052
## 22499 -122.4438 37.72965
## 22531 -122.4141 37.77608
## 22533 -122.4409 37.74887
## 22549 -122.4431 37.71275
## 22554 -122.4009 37.76461
## 22600 -122.4362 37.74701
## 22607 -122.4753 37.71516
## 22673 -122.4131 37.74787
## 22739 -122.3919 37.78868
## 22801 -122.4709 37.75504
## 22802 -122.4014 37.77974
## 22820 -122.4887 37.77819
## 22870 -122.4367 37.76372
## 22871 -122.4507 37.72158
## 22896 -122.4342 37.74955
## 22924 -122.4460 37.71426
## 22931 -122.4210 37.74497
## 22951 -122.4948 37.77410
## 23000 -122.4254 37.76334
## 23036 -122.4163 37.80401
## 23045 -122.3996 37.77438
## 23056 -122.4753 37.71487
## 23074 -122.3696 37.72900
## 23117 -122.4223 37.79267
## 23135 -122.3942 37.77181
## 23154 -122.4648 37.75510
## 23172 -122.4415 37.78044
## 23173 -122.4835 37.78257
## 23251 -122.4119 37.71432
## 23267 -122.4263 37.71022
## 23308 -122.4769 37.76468
## 23356 -122.4613 37.71719
## 23400 -122.4689 37.77570
## 23433 -122.4753 37.71553
## 23467 -122.3928 37.78205
## 23478 -122.4855 37.75523
## 23486 -122.4276 37.72233
## 23503 -122.4108 37.76534
## 23519 -122.4454 37.76390
## 23535 -122.4639 37.76393
## 23540 -122.4599 37.72230
## 23542 -122.4860 37.78680
## 23598 -122.3910 37.78938
## 23601 -122.3922 37.78607
## 23721 -122.4452 37.71600
## 23782 -122.4012 37.75378
## 23790 -122.5038 37.75594
## 23828 -122.3957 37.73636
## 23881 -122.3906 37.75501
## 23895 -122.4517 37.72883
## 23896 -122.4011 37.77047
## 23913 -122.4734 37.75749
## 23918 -122.4676 37.71048
## 23919 -122.4936 37.77938
## 23991 -122.4418 37.73477
## 24130 -122.4544 37.76267
## 24185 -122.4891 37.77817
## 24223 -122.4194 37.75596
## 24282 -122.4934 37.73923
## 24314 -122.4577 37.74434
## 24326 -122.3880 37.73296
## 24374 -122.3696 37.72900
## 24393 -122.4372 37.75934
## 24401 -122.4517 37.73329
## 24402 -122.4711 37.74663
## 24504 -122.4536 37.70915
## 24546 -122.4544 37.78545
## 24554 -122.4395 37.75941
## 24574 -122.3997 37.75216
## 24607 -122.4836 37.77285
## 24639 -122.4824 37.75738
## 24659 -122.3920 37.72865
## 24674 -122.4760 37.78306
## 24700 -122.4138 37.73519
## 24716 -122.4922 37.77477
## 24724 -122.4374 37.75893
## 24742 -122.4704 37.78430
## 24781 -122.4137 37.79099
## 24794 -122.4016 37.75125
## 24816 -122.4009 37.72799
## 24840 -122.4704 37.74936
## 24872 -122.4236 37.77800
## 24896 -122.4441 37.73437
## 24922 -122.4446 37.77193
## 24940 -122.4192 37.72509
## 24949 -122.4356 37.72596
## 24966 -122.3919 37.78868
## 24983 -122.4458 37.73512
## 24986 -122.4668 37.77558

Compare the SPDF to DF

str(sfhomes15) # the data frame

str(sfhomes15_sp) # the data frame

Compare the SPDF to DF

str(sfhomes15) # the data frame

## 'data.frame':    835 obs. of  19 variables:
##  $ FiscalYear      : int  2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##  $ SalesDate       : chr  "2015-08-21" "2015-08-13" "2015-12-29" "2015-08-26" ...
##  $ Address         : chr  "0000 2760 19TH                AV0015" "0000 0560AMISSOURI            ST0000" "0000 0718 LONG BRIDGE         ST1202" "0000 0103 SUMMIT              WY0000" ...
##  $ YearBuilt       : int  1979 2003 2016 NA 2015 1961 1973 1900 NA 2015 ...
##  $ NumBedrooms     : int  2 2 2 3 3 3 3 2 0 2 ...
##  $ NumBathrooms    : int  2 2 2 3 3 3 3 2 0 2 ...
##  $ NumRooms        : int  5 5 5 6 7 7 7 5 0 0 ...
##  $ NumStories      : int  0 1 1 0 1 14 2 1 0 1 ...
##  $ NumUnits        : int  1 1 1 1 1 1 1 1 0 1 ...
##  $ AreaSquareFeet  : int  1595 1191 1346 2133 1266 1840 2075 1256 0 1520 ...
##  $ ImprovementValue: int  432500 701280 1512093 582108 850000 1154846 484544 810169 486872 962500 ...
##  $ LandValue       : int  432500 701280 748900 873161 850000 1154846 1130604 1890395 1136036 962500 ...
##  $ Neighborhood    : chr  "West of Twin Peaks" "Potrero Hill" "Mission Bay" "Lakeshore" ...
##  $ Location        : chr  "(37.7360097396496, -122.474067310226)" "(37.759197817252, -122.396516184449)" "(37.7718055392903, -122.394186974738)" "(37.7151803909173, -122.475425717555)" ...
##  $ SupeDistrict    : int  7 10 6 7 9 3 8 2 7 6 ...
##  $ totvalue        : int  865000 1402560 2260993 1455269 1700000 2309692 1615148 2700564 1622908 1925000 ...
##  $ SalesYear       : int  2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##  $ lat             : num  37.7 37.8 37.8 37.7 37.8 ...
##  $ lon             : num  -122 -122 -122 -122 -122 ...

SPDF

str(sfhomes15_sp) # the SPDF

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  835 obs. of  17 variables:
##   .. ..$ FiscalYear      : int [1:835] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##   .. ..$ SalesDate       : chr [1:835] "2015-08-21" "2015-08-13" "2015-12-29" "2015-08-26" ...
##   .. ..$ Address         : chr [1:835] "0000 2760 19TH                AV0015" "0000 0560AMISSOURI            ST0000" "0000 0718 LONG BRIDGE         ST1202" "0000 0103 SUMMIT              WY0000" ...
##   .. ..$ YearBuilt       : int [1:835] 1979 2003 2016 NA 2015 1961 1973 1900 NA 2015 ...
##   .. ..$ NumBedrooms     : int [1:835] 2 2 2 3 3 3 3 2 0 2 ...
##   .. ..$ NumBathrooms    : int [1:835] 2 2 2 3 3 3 3 2 0 2 ...
##   .. ..$ NumRooms        : int [1:835] 5 5 5 6 7 7 7 5 0 0 ...
##   .. ..$ NumStories      : int [1:835] 0 1 1 0 1 14 2 1 0 1 ...
##   .. ..$ NumUnits        : int [1:835] 1 1 1 1 1 1 1 1 0 1 ...
##   .. ..$ AreaSquareFeet  : int [1:835] 1595 1191 1346 2133 1266 1840 2075 1256 0 1520 ...
##   .. ..$ ImprovementValue: int [1:835] 432500 701280 1512093 582108 850000 1154846 484544 810169 486872 962500 ...
##   .. ..$ LandValue       : int [1:835] 432500 701280 748900 873161 850000 1154846 1130604 1890395 1136036 962500 ...
##   .. ..$ Neighborhood    : chr [1:835] "West of Twin Peaks" "Potrero Hill" "Mission Bay" "Lakeshore" ...
##   .. ..$ Location        : chr [1:835] "(37.7360097396496, -122.474067310226)" "(37.759197817252, -122.396516184449)" "(37.7718055392903, -122.394186974738)" "(37.7151803909173, -122.475425717555)" ...
##   .. ..$ SupeDistrict    : int [1:835] 7 10 6 7 9 3 8 2 7 6 ...
##   .. ..$ totvalue        : int [1:835] 865000 1402560 2260993 1455269 1700000 2309692 1615148 2700564 1622908 1925000 ...
##   .. ..$ SalesYear       : int [1:835] 2015 2015 2015 2015 2015 2015 2015 2015 2015 2015 ...
##   ..@ coords.nrs : int [1:2] 19 18
##   ..@ coords     : num [1:835, 1:2] -122 -122 -122 -122 -122 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:835] "24" "35" "76" "79" ...
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   ..@ bbox       : num [1:2, 1:2] -122.5 37.7 -122.4 37.8
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA

SPDF Objects

You can see from str(sfhomes) that a SPDF object is a collection of slots or components. The key ones are:

@data data frame of attributes that describe each location
@coords the coordinates for each geometric object - here points
@bbox the min and max bounding coordinates
@proj4string the coordinate reference system defintion as a string

Explore the object in the Environment window

SPDF Slots

Review the output of each of these:

summary(sfhomes15_sp)
head(sfhomes15_sp@data)
class(sfhomes15_sp@data)

sfhomes15_sp@bbox
bbox(sfhomes15_sp)

head(sfhomes15_sp@coords)
head(sfhomes15_sp$lat)
head(sfhomes15_sp$lon)

sfhomes15_sp@proj4string
proj4string(sfhomes15_sp)

Take a closer look

Look at sfhomes15_sp in the environment window

What’s missing

Are all the columns that were present in sfhomes15 also in sfhomes15_sp?

Is there a slot in sfhomes15_sp without data?

What is the CRS of the data?

proj4string(sfhomes15_sp) # get a CRS object

## [1] NA

Map a `SpatialPointsDataFrame`

plot(sfhomes15_sp)  # using sp::plot

So…..?

Recap

We created great maps of sfhomes point data with ggplot and ggmap.

Then we created a simple map of the SPDF with plot.

We aren’t seeing the value of sp objects just yet.

Let’s add more geospatial data

`rgdal`

rgdal is an R port of the powerful and widely used GDAL library.

It is the most commonly used R library for importing and exporting spatial data.

OGR: for vector data: readOGR() and writeOGR()
GDAL for raster data: readGDAL() and writeGDAL()

`rgdal`

library(rgdal)

## Warning: package 'rgdal' was built under R version 3.4.3

## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.4/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-5

# See what file types are supported by rgdal drivers
# ogrDrivers()$name

Getting help

gdal.org

?readOGR

For more info on working with rgdal to load different types of spatial data in R see this excellent tutorial by Zev Ross.

Read in a Shapefile with the boundary of San Francisco

Take a look at the file(s)

dir("data", pattern="sf_boundary")

## [1] "sf_boundary.dbf" "sf_boundary.prj" "sf_boundary.shp" "sf_boundary.shx"

Read in Shapefile

sfboundary <- readOGR(dsn="data",layer="sf_boundary")

## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "sf_boundary"
## with 1 features
## It has 3 fields

# or
# sfboundary <- readOGR("data","sf_boundary")
# but not
#sfboundary <- readOGR(dsn="data/",layer="sf_boundary")

Check out the data structure

What type of sp object is sfboundary? How many features are in the object?

class(sfboundary)
str(sfboundary) 
head(sfboundary@data)

Explore the object in the Envi window

Make a quick plot of `sfboundary`

How?

Make a quick plot of `sfboundary`

plot(sfboundary)

Map `sfboundary` & `sfhomes15_sp`

plot(sfboundary)
points(sfhomes15_sp, col="red")

Map `sfboundary` & `sfhomes15_sp`

Where are the points? What’s wrong?

plot(sfboundary)
points(sfhomes15_sp, col="red")

What’s Wrong?

Compare the CRSs, are they the same?

proj4string(sfboundary)
proj4string(sfhomes15_sp)

Compare the CRSs, are they the same?

proj4string(sfboundary)

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

proj4string(sfhomes15_sp)

## [1] NA

Compare the coordinate values

sfboundary@bbox

##         min       max
## x  542696.6  556659.9
## y 4173563.7 4185088.6

sfhomes15_sp@bbox

##            min        max
## lon -122.51072 -122.36964
## lat   37.70854   37.80612

CRS Problems

The #1 reason…

CRS Definitions

All sp objects should have a defined CRS

If not, one needs to be assigned to the object

This is called defining a projection
This doesn’t change the coordinates

CRS Transformations

All sp objects should have the same CRS.

This CRS should be defined in the proj4string slot.

When they don’t, they need to be transformed to a common CRS.
This is called a projection transformation,
- or projecting or reprojection.

Projection transformation returns a new spatial object with the transformed coordinates

CRS Definitions and Transformations

GIS software - or software for working with geospatial data - will have a database of definitions for thousands of Earth referenced coordinate systems

and methods for using those CRS definitions to define or transform a CRS.

Defining a CRS

So, do we need to define the CRS of either

sfboundary or
sfhomes15_sp?

If yes, which one? How can you tell?

Defining a CRS

We need to know

what the data’s CRS is
how to define it for the sp object

What is the CRS

bbox(sfhomes15_sp)

##            min        max
## lon -122.51072 -122.36964
## lat   37.70854   37.80612

Defining a CRS

Use the sp function CRS() to create a CRS object

Then, use the sp::proj4string() to assign it

# Define the CRS for sfhomes15_sp as WGS84

proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
                               +ellps=WGS84 +datum=WGS84 +no_defs")

Referencing the WGS84 CRS

The WGS84 proj4 string

CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

The WGS84 epsg code

CRS("+init=epsg:4326")

EPSG refers to the group that defined these

European Petroleum Survey Group

EPSG Example

EPSG codes are more commonly used. Can you guess why?

# use an EPSG code for WGS84
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326") 

# or enter the parameter string
# proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
#                               +ellps=WGS84 +datum=WGS84 +no_defs")

Proj4 CRS Definitions

Proj4 is the standard library for defining and transforming map projections and coordinate reference systems.

Here is an example file of proj4 CRS definitions

Check it

Once you set the CRS you can check it with the same function.

proj4string(sfhomes15_sp)

## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

Recap

proj4string() can get or set the CRS

CRS() is used to define the CRS of an sp object

The following are equivalent

# Define the CRS with the parameter string
proj4string(sfhomes15_sp) <- CRS("+proj=longlat 
                              +ellps=WGS84 +datum=WGS84 +no_defs")

# Define the CRS with an EPSG code for WGS84
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326")

epsg codes are used as short-hand for CRSs definitions

Compare the CRSs, again

Are they both defined? Are they the same?

proj4string(sfboundary)

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

proj4string(sfhomes15_sp)

## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"

proj4string(sfboundary) == proj4string(sfhomes15_sp)

## [1] FALSE

Transforming a CRS

Transform the CRS

Use sp function spTransform

Requires as input:

a sp object to transform with a defined CRS
a target CRS

Outputs a new spatial object with coordinate data in the target CRS

It does not change the input data

Transform `sfboundary`

Let’s transform the CRS of sfboundary to be the same as that of sfhomes15_sp.

sfboundary_lonlat <- spTransform(sfboundary, CRS("+init=epsg:4326"))

Why not the other way around?

Did it work?

Do they match now?

How will we know?

Do the CRSs match?

proj4string(sfhomes15_sp) == proj4string(sfboundary_lonlat)

## [1] FALSE

Overlay the data in space

plot(sfboundary_lonlat)
points(sfhomes15_sp, col="red")

Overlay the data in space

Woo-hoo!

So….

We can define the CRS of sp objects and transform them to the same CRS for mapping and spatial analysis.

Recap

# Define the CRS
proj4string(sfhomes15_sp) <- CRS("+init=epsg:4326")

#Transform the CRS
sfboundary_lonlat <- spTransform(sfboundary, CRS("+init=epsg:4326"))

Save `sfboundary_lonlat`

Use writeOGR to save sfboundary_lonlat to a new shapefile

See ?writeOGR for help

# write transformed data to a new shapefile 
writeOGR(sfboundary_lonlat, 
          dsn = "data", 
          layer = "sfbounary_lonlat", 
          driver="ESRI Shapefile")

# is it there?
dir("data")

Projections, CRS, oh my!

We want all data in the same CRS

Which one is best?

Common CRS Codes

Geographic CRSs

4326 Geographic, WGS84 (default for lon/lat)
4269 Geographic, NAD83 (USA Fed agencies like Census)

Projected CRSs

5070 USA Contiguous Albers Equal Area Conic
3310 CA ALbers Equal Area
26910 UTM Zone 10, NAD83 (Northern Cal)
3857 Web Mercator (web maps)

Finding CRS Codes

See http://spatialreference.org/

Use this site to find EPSG codes and proj4 CRS strings

Challenge

Use coordinates to convert the landmarks data frame to a SpatialPointsDataFrame
Use proj4string to define it’s CRS as 3857 Web Mercator
Use spTransform to transform the landmarks data to geographic coords (WGS84)
Add the landmarks to the last map we created.

Challenge - solution

# Convert the DF to a SPDF
coordinates(landmarks) <-c("X","Y")
# Define the CRS with an EPSG code for Web mercator
proj4string(landmarks) <- CRS("+init=epsg:3857") 
# Transform the CRS to WGS84
landmarks_lonlat <- spTransform(landmarks, CRS("+init=epsg:4326"))
# map it
plot(sfboundary_lonlat)
points(sfhomes15_sp, col="red")
points(landmarks_lonlat, col="green")

Challenge - solution

Challenge

Use readOGR to load the shapefile “data/sf_highways.shp”
Check its data type with class()
Check its CRS with proj4string()
Use spTransform if need to project it to WGS84
Recreate the last map and add the highways as black lines

Challenge - solution

highways <- readOGR(dsn="data", layer="sf_highways")
class(highways)
proj4string(highways)
highways_lonlat <- spTransform(highways, CRS("+init=epsg:4326"))

plot(sfboundary_lonlat)
lines(highways_lonlat, col="black")
points(sfhomes15_sp, col="red")
points(landmarks_lonlat, col="green")

Challenge - solution

## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "sf_highways"
## with 246 features
## It has 5 fields

## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"

## [1] "+proj=utm +zone=10 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

# QUESTIONS?

Recap

Geospatial Data in R
CSV > Data Frame > ggplot/ggmap
sp library - spatial objects and methods in R
GDAL libary - readOGR, writeOGR
Convert Data Frame to SPDF
CRS definitions and transformations
Map Overlays

Geospatial Data in R, pt1

Patty Frontiera

Spring 2018

Workshop Prep

Geospatial Data in R

Workshop Goals

About Me

About you

Follow Along

Geographic Data

Geographic Data

Geospatial data

Coordinate Reference Systems

Coordinate Reference Systems

WGS84

Geospatial data are powerful!

Spatial Data

Types of Spatial Data

Types of Spatial Data

Vector Data

Raster Data

Softare for working with Geospatial Data

Geospatial data require

GIS

Types of GIS Software

Why R for Geospatial Data?

Why R for Geospatial Data?

Geospatial Data in R

Prep reminder

Geospatial Data in R

About the Sample Data

Load the CSV file into a data frame

Explore the data

Questions:

Plot of points

Plot of points

ggplot2

ggplot2

ggplot2

Maps with ggplot2

Maps with ggplot2

Coord_map Option

coord_map option

Map Projections

Many Map Projections & Projected CRSs

Different Projected CRSs

Map points symbolized by totvalue

Map points symbolized by totvalue

Data Order

Data Order

Challenge

Decreasing order by totvalue

More ggplot Goodness

More ggplot Goodness

ggmap

ggmap extends ggplot

ggmap

ggmap

Register Google Maps Key

Re-install ggmap

get_map

get_map

Display the ggmap of SF

Display the ggmap of SF

ggmap with point overlay

ggmap with point overlay

Customize get_map

Customize get_map Extent

Challenge

Map the points with sf_basemap

Map Overlays

Bart Stations and Landmarks

Add BART stations to map

Add BART stations to map

Geocoding with GGMAP

Geocoding

Geocode SF Police Stations

Geocode the Police Stations

Single column address

Geocode

`ggplot2`

`ggplot2`

`ggplot2`

Maps with `ggplot2`

Maps with `ggplot2`

Map points symbolized by `totvalue`

Map points symbolized by `totvalue`

Decreasing order by `totvalue`

More `ggplot` Goodness

More `ggplot` Goodness

`ggmap` extends `ggplot`

Re-install `ggmap`

`ggmap` with point overlay

`ggmap` with point overlay

Customize `get_map`

Customize `get_map` Extent

Map the points with `sf_basemap`

The `sp` Package

`sp` package

`sp::coordinates()`

Map a `SpatialPointsDataFrame`

`rgdal`

`rgdal`

Make a quick plot of `sfboundary`

Make a quick plot of `sfboundary`

Map `sfboundary` & `sfhomes15_sp`

Map `sfboundary` & `sfhomes15_sp`